Summary:

Two avenues of analysis:

  1. First we see how \(\Delta\)load affects \(\Delta\)ArcLength^2 and \(\Delta\)frequency^2, accounting for bee size and treatment order
  2. We check to see how \(\Delta\)frequency and \(\Delta\)ArcLength are associated with \(\Delta\)MetabolicRate, while accounting for \(\Delta\)load

Measured variables:

Calculated variables:

Install packages and read data

ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS",'reshape2', 
              'influence.ME', 'sjPlot', "effects", 'visreg', 'viridis')
ipak(packages)
##      ggplot2          car         lme4       gsheet         MASS 
##         TRUE         TRUE         TRUE         TRUE         TRUE 
##     reshape2 influence.ME       sjPlot      effects       visreg 
##         TRUE         TRUE         TRUE         TRUE         TRUE 
##      viridis 
##         TRUE
# read in data -- google sheet called "Final Bee Resp Data"
url = 'https://docs.google.com/spreadsheets/d/1wT-QxSJJElhhJcIXlg2hDpFKHbLyFGuqNNj2iYvA8Vo/edit?usp=sharing'
bdta <- data.frame(gsheet2tbl(url))
summary(bdta)
##     BeeID               order      Treatment            Mstarved     
##  Length:60          Min.   :1.0   Length:60          Min.   :0.0767  
##  Class :character   1st Qu.:1.0   Class :character   1st Qu.:0.1213  
##  Mode  :character   Median :1.5   Mode  :character   Median :0.1378  
##                     Mean   :1.5                      Mean   :0.1411  
##                     3rd Qu.:2.0                      3rd Qu.:0.1670  
##                     Max.   :2.0                      Max.   :0.1909  
##        M2               MF             Itspan               S            
##  Min.   :0.1062   Min.   :0.0969   Min.   :0.003860   Min.   :4.010e-05  
##  1st Qu.:0.1808   1st Qu.:0.1762   1st Qu.:0.004400   1st Qu.:5.590e-05  
##  Median :0.2191   Median :0.2107   Median :0.004620   Median :6.130e-05  
##  Mean   :0.2245   Mean   :0.2154   Mean   :0.004594   Mean   :6.141e-05  
##  3rd Qu.:0.2646   3rd Qu.:0.2576   3rd Qu.:0.004870   3rd Qu.:6.750e-05  
##  Max.   :0.3550   Max.   :0.3440   Max.   :0.005180   Max.   :7.810e-05  
##       MetR             freq            amp               L           
##  Min.   : 5.507   Min.   :163.9   Min.   : 96.88   Min.   :0.008691  
##  1st Qu.: 8.935   1st Qu.:173.0   1st Qu.:114.16   1st Qu.:0.010570  
##  Median :10.666   Median :179.1   Median :125.96   Median :0.011017  
##  Mean   :10.600   Mean   :178.9   Mean   :123.31   Mean   :0.010911  
##  3rd Qu.:12.297   3rd Qu.:185.3   3rd Qu.:133.61   3rd Qu.:0.011444  
##  Max.   :16.475   Max.   :194.7   Max.   :144.78   Max.   :0.012118

Calculate variables

bdta <- within(bdta, {
     MT = (M2 + MF)/2
     load = MT - Mstarved
     perLoad = load / Mstarved * 100
     arcL = (0.75 * L) * (amp * (pi / 180))
     U = arcL * freq * 2
     freq2 = freq^2
     arcL2 = arcL^2
     frce = U^2 * S
     
     # convert to factor variables
     order = as.factor(as.character(order))
     Treatment = as.factor(as.character(Treatment))
     BeeID = as.factor(as.character(BeeID))
})

Create a new dataframe that calculates the changes for each bee

# convert from long to wide
newDF <- data.frame()
colsTocalc = c("order", "M2", "MF", "MetR", "freq", "amp", "load", "MT", "perLoad", "frce", "arcL2", "freq2", "U", "arcL")
for(varb in colsTocalc){
     data_wide <- dcast(bdta, BeeID + Mstarved + S + Itspan + L ~ 
                             Treatment, value.var=c(varb))
     colnames(data_wide)[6:7] = paste(varb, colnames(data_wide)[6:7], sep = "_")
     if(varb == colsTocalc[1]){
       newDF <- data_wide   
     }
     else newDF <- merge(newDF, data_wide, all.y = TRUE)    
}

head(newDF)
##   BeeID Mstarved        S  Itspan          L order_H order_L   M2_H   M2_L
## 1   E01   0.1577 6.46e-05 0.00464 0.01116843       2       1 0.2634 0.1895
## 2   E03   0.1732 6.92e-05 0.00487 0.01144364       1       2 0.3106 0.2105
## 3   E05   0.1755 7.06e-05 0.00497 0.01168730       2       1 0.2776 0.2006
## 4   E09   0.1135 4.88e-05 0.00440 0.00986884       1       2 0.2282 0.1425
## 5   E10   0.1269 5.91e-05 0.00462 0.01060036       1       2 0.2166 0.1460
## 6   E11   0.1350 5.88e-05 0.00464 0.01070018       1       2 0.2223 0.1599
##     MF_H   MF_L   MetR_H    MetR_L   freq_H   freq_L    amp_H     amp_L
## 1 0.2559 0.1812 12.54241 10.393483 186.9973 179.3308 128.7145 111.49875
## 2 0.3079 0.2023 11.78387  9.214622 174.9374 167.8412 144.7784 119.98953
## 3 0.2742 0.1975 12.20102  8.981412 178.8896 166.2385 125.6778 116.22053
## 4 0.2084 0.1390  9.47114  7.581104 191.8348 190.3535 139.3246 103.30703
## 5 0.2099 0.1422 11.03371  7.159840 186.1762 169.5451 128.8257  96.88003
## 6 0.2179 0.1579 10.23580  6.806696 183.0259 168.2223 128.7916 116.18262
##    load_H  load_L    MT_H    MT_L perLoad_H perLoad_L      frce_H
## 1 0.10195 0.02765 0.25965 0.18535  64.64807  17.53329 0.003199482
## 2 0.13605 0.03320 0.30925 0.20640  78.55081  19.16859 0.003984230
## 3 0.10040 0.02355 0.27590 0.19905  57.20798  13.41880 0.003340857
## 4 0.10480 0.02725 0.21830 0.14075  92.33480  24.00881 0.002327019
## 5 0.08635 0.01720 0.21325 0.14410  68.04571  13.55398 0.002618299
## 6 0.08510 0.02390 0.22010 0.15890  63.03704  17.70370 0.002563877
##        frce_L      arcL2_H      arcL2_L  freq2_H  freq2_L      U_H
## 1 0.002208025 0.0003540922 0.0002657061 34967.99 32159.55 7.037583
## 2 0.002519159 0.0004703414 0.0003230668 30603.08 28170.67 7.587858
## 3 0.002467172 0.0003696773 0.0003161342 32001.50 27635.24 6.879020
## 4 0.001259710 0.0003239404 0.0001781022 36800.61 36234.44 6.905419
## 5 0.001228018 0.0003195386 0.0001807120 34661.58 28745.53 6.656039
## 6 0.001762569 0.0003254128 0.0002648146 33498.49 28298.73 6.603283
##        U_L     arcL_H     arcL_L
## 1 5.846363 0.01881734 0.01630050
## 2 6.033576 0.02168736 0.01797406
## 3 5.911496 0.01922699 0.01778016
## 4 5.080722 0.01799834 0.01334549
## 5 4.558361 0.01787564 0.01344292
## 6 5.475004 0.01803920 0.01627313

Calculate \(\Delta\) variables

newDF <- within(newDF, {
     deltaPercLoad = perLoad_H - perLoad_L
     avgPercLoad = (perLoad_H + perLoad_L) / 2
     deltaMetR = MetR_H - MetR_L
     deltaFrq2 = freq2_H - freq2_L
     deltaArcL = arcL_H - arcL_L
     deltaArcL2 = arcL2_H - arcL2_L
     deltaFreq2Perc = deltaFrq2 / deltaPercLoad
     deltaLoad = scale(load_H - load_L, center = TRUE, scale = FALSE)
     dLoad_nonCent = load_H - load_L
     deltaLoad2 <- deltaLoad^2
})


plot(newDF$deltaArcL2 ~ newDF$deltaFrq2)

# compare with Susie's calculations



b2 <- read.csv("~/Desktop/b2.csv")
b2 <- b2[!is.na(b2$deltaperload), ]
b2 <- b2[order(b2$BeeID, b2$Treatment), ]

plot(b2$AverageperLoad, newDF$avgPercLoad)
abline(0,1)

plot(b2$deltaload, newDF$deltaLoad)
abline(0,1)

Use PCA to combine bee sizes into a single predictor

# principle components
aa = prcomp(newDF[, c("Mstarved", "Itspan", "S", "L")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size
## Importance of components:
##                           PC1     PC2     PC3     PC4
## Standard deviation     1.9422 0.37234 0.22746 0.19332
## Proportion of Variance 0.9431 0.03466 0.01293 0.00934
## Cumulative Proportion  0.9431 0.97772 0.99066 1.00000
biplot(aa) # shows that all three size measurement are correlated

# note, I changed the signs of the predictions so that higher PC1 values 
# correspond to bigger bees
p1 = -predict(aa)[,1] 

# add PC1 scores to dataset
newDF$size_pc1 = p1
newDF$size_pc1_2 = newDF$size_pc1^2

# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(newDF[, c("Mstarved", "Itspan", "S",  "L",  "size_pc1")])

1. First we see how \(\Delta\)load affects \(\Delta\)ArcLength accounting for bee size and treatment order

This is the model selection procedure that was used: 1. Fit a large model with all two-way interactions and squared terms 2. Remove non-significant predictors, starting with interactions, squared terms, and then main effects

# reformat order so that it is more interpretable
library(plyr)
newDF$order_1 <- mapvalues(newDF$order_H, from = c(2, 1), to = c("loadedSecond", "loadedFirst"))

# fit full model
m1 <- lm(deltaArcL2 ~  (deltaLoad +  size_pc1 + order_1)^2  + 
              size_pc1_2  + deltaLoad2, data = newDF)
summary(m1)
## 
## Call:
## lm(formula = deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 + 
##     size_pc1_2 + deltaLoad2, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.326e-05 -1.951e-05  3.244e-06  1.635e-05  4.401e-05 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    9.818e-05  1.300e-05   7.551 2.05e-07 ***
## deltaLoad                      1.859e-03  1.178e-03   1.578    0.130    
## size_pc1                      -1.804e-05  1.263e-05  -1.428    0.168    
## order_1loadedSecond            7.310e-07  1.978e-05   0.037    0.971    
## size_pc1_2                    -3.660e-06  4.905e-06  -0.746    0.464    
## deltaLoad2                     3.807e-02  4.366e-02   0.872    0.393    
## deltaLoad:size_pc1            -2.300e-06  8.387e-04  -0.003    0.998    
## deltaLoad:order_1loadedSecond -6.487e-04  2.225e-03  -0.292    0.774    
## size_pc1:order_1loadedSecond   1.204e-05  2.099e-05   0.574    0.572    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.03e-05 on 21 degrees of freedom
## Multiple R-squared:  0.575,  Adjusted R-squared:  0.413 
## F-statistic: 3.551 on 8 and 21 DF,  p-value: 0.009387
m2 <- update(m1, .~. - deltaLoad:size_pc1)
anova(m1, m2)
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 + 
##     deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1 + size_pc1:order_1
##   Res.Df        RSS Df   Sum of Sq  F Pr(>F)
## 1     21 1.9281e-08                         
## 2     22 1.9281e-08 -1 -6.9071e-15  0 0.9978
summary(m2)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + deltaLoad:order_1 + size_pc1:order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.327e-05 -1.953e-05  3.263e-06  1.634e-05  4.402e-05 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    9.817e-05  1.252e-05   7.840 8.24e-08 ***
## deltaLoad                      1.860e-03  9.935e-04   1.872   0.0745 .  
## size_pc1                      -1.806e-05  7.886e-06  -2.290   0.0320 *  
## order_1loadedSecond            7.471e-07  1.846e-05   0.040   0.9681    
## size_pc1_2                    -3.672e-06  1.745e-06  -2.104   0.0470 *  
## deltaLoad2                     3.800e-02  3.336e-02   1.139   0.2669    
## deltaLoad:order_1loadedSecond -6.509e-04  2.022e-03  -0.322   0.7505    
## size_pc1:order_1loadedSecond   1.209e-05  1.187e-05   1.019   0.3194    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.96e-05 on 22 degrees of freedom
## Multiple R-squared:  0.575,  Adjusted R-squared:  0.4397 
## F-statistic: 4.251 on 7 and 22 DF,  p-value: 0.004168
m3 <- update(m2, .~. - deltaLoad:order_1)

anova(m2, m3)
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1 + size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     size_pc1:order_1
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     22 1.9281e-08                             
## 2     23 1.9371e-08 -1 -9.0841e-11 0.1037 0.7505
summary(m3)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + size_pc1:order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.414e-05 -1.856e-05  3.421e-06  1.715e-05  4.272e-05 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   9.768e-05  1.218e-05   8.017 4.13e-08 ***
## deltaLoad                     1.625e-03  6.585e-04   2.467   0.0215 *  
## size_pc1                     -1.728e-05  7.355e-06  -2.349   0.0278 *  
## order_1loadedSecond           3.683e-06  1.573e-05   0.234   0.8170    
## size_pc1_2                   -3.816e-06  1.654e-06  -2.308   0.0304 *  
## deltaLoad2                    4.504e-02  2.467e-02   1.825   0.0809 .  
## size_pc1:order_1loadedSecond  1.003e-05  9.807e-06   1.023   0.3169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.902e-05 on 23 degrees of freedom
## Multiple R-squared:  0.573,  Adjusted R-squared:  0.4615 
## F-statistic: 5.143 on 6 and 23 DF,  p-value: 0.001756
m4 <- update(m3, .~. - size_pc1:order_1)
anova(m3, m4) 
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     23 1.9371e-08                             
## 2     24 2.0253e-08 -1 -8.8175e-10 1.0469 0.3169
summary(m4)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.427e-05 -1.874e-05 -5.000e-09  1.783e-05  4.254e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.024e-04  1.128e-05   9.079 3.13e-09 ***
## deltaLoad            1.639e-03  6.590e-04   2.486   0.0203 *  
## size_pc1            -1.180e-05  5.042e-06  -2.340   0.0279 *  
## order_1loadedSecond -5.411e-07  1.519e-05  -0.036   0.9719    
## size_pc1_2          -2.627e-06  1.178e-06  -2.231   0.0353 *  
## deltaLoad2           2.613e-02  1.636e-02   1.597   0.1233    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.905e-05 on 24 degrees of freedom
## Multiple R-squared:  0.5535, Adjusted R-squared:  0.4605 
## F-statistic: 5.951 on 5 and 24 DF,  p-value: 0.001027
m5 <- update(m4, .~. - deltaLoad2)
anova(m4,m5) 
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     24 2.0253e-08                             
## 2     25 2.2406e-08 -1 -2.1529e-09 2.5511 0.1233
summary(m5)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2, 
##     data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.782e-05 -1.748e-05  9.200e-07  1.678e-05  4.278e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.074e-04  1.117e-05   9.623 6.93e-10 ***
## deltaLoad            2.095e-03  6.122e-04   3.422  0.00215 ** 
## size_pc1            -1.388e-05  5.018e-06  -2.767  0.01050 *  
## order_1loadedSecond  2.783e-07  1.565e-05   0.018  0.98595    
## size_pc1_2          -2.037e-06  1.152e-06  -1.768  0.08930 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.994e-05 on 25 degrees of freedom
## Multiple R-squared:  0.5061, Adjusted R-squared:  0.427 
## F-statistic: 6.403 on 4 and 25 DF,  p-value: 0.001086
m6 <- update(m5, .~.  - size_pc1_2)
anova(m6, m5) 
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)  
## 1     26 2.5207e-08                              
## 2     25 2.2406e-08  1 2.8009e-09 3.1251 0.0893 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m6)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.154e-05 -1.781e-05 -4.100e-07  1.849e-05  4.983e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          9.844e-05  1.034e-05   9.525 5.79e-10 ***
## deltaLoad            2.148e-03  6.359e-04   3.378  0.00231 ** 
## size_pc1            -1.321e-05  5.205e-06  -2.539  0.01745 *  
## order_1loadedSecond  3.224e-06  1.618e-05   0.199  0.84364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.114e-05 on 26 degrees of freedom
## Multiple R-squared:  0.4443, Adjusted R-squared:  0.3802 
## F-statistic: 6.929 on 3 and 26 DF,  p-value: 0.001402
m7 <- update(m6, .~. - order_1)
anova(m7, m6)
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)
## 1     27 2.5245e-08                            
## 2     26 2.5207e-08  1 3.8477e-11 0.0397 0.8436
summary(m7)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.002e-04  5.583e-06  17.940  < 2e-16 ***
## deltaLoad    2.059e-03  4.439e-04   4.638 8.06e-05 ***
## size_pc1    -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666
anova(m7, update(m7, .~. - deltaLoad)) # p-value for delta load
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ size_pc1
##   Res.Df        RSS Df   Sum of Sq      F    Pr(>F)    
## 1     27 2.5245e-08                                    
## 2     28 4.5359e-08 -1 -2.0114e-08 21.512 8.055e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m7, update(m7, .~. - size_pc1)) # p-value for size
## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad
##   Res.Df        RSS Df   Sum of Sq      F   Pr(>F)   
## 1     27 2.5245e-08                                  
## 2     28 3.4619e-08 -1 -9.3737e-09 10.025 0.003808 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m7)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.002e-04  5.583e-06  17.940  < 2e-16 ***
## deltaLoad    2.059e-03  4.439e-04   4.638 8.06e-05 ***
## size_pc1    -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666
# refit model with non-centered version of deltaLoad
# this model will have a different intercept, but same p-values for slopes
summary(update(m7, .~. - deltaLoad + dLoad_nonCent)) # final model for paper
## 
## Call:
## lm(formula = deltaArcL2 ~ size_pc1 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -5.822e-05  3.460e-05  -1.683  0.10397    
## size_pc1      -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## dLoad_nonCent  2.059e-03  4.439e-04   4.638 8.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666

model diagnostics

par(mfrow = c(2,2))
plot(m7, which = 1:4) # no glaring violations, though there are a few fairly influential observations

par(mfrow = c(1,1))

car::vif(m7) 
## deltaLoad  size_pc1 
##  1.841066  1.841066

Model visualization

summary(m7)
## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.002e-04  5.583e-06  17.940  < 2e-16 ***
## deltaLoad    2.059e-03  4.439e-04   4.638 8.06e-05 ***
## size_pc1    -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666
# calculate partial residuals for deltaLoad
# these are the  residuals, minus the effect of detlaLoad
y <- residuals(m7, type = 'partial')[, "deltaLoad"]

# plot partial residuals with base R plotting
plot(x = newDF$deltaLoad, y = y)

# double check to make sure the slope for partial residual plots are the 
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent)) 
## 
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1.584e-04  2.531e-05  -6.257 9.20e-07 ***
## newDF$dLoad_nonCent  2.059e-03  3.213e-04   6.409 6.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.003e-05 on 28 degrees of freedom
## Multiple R-squared:  0.5946, Adjusted R-squared:  0.5801 
## F-statistic: 41.07 on 1 and 28 DF,  p-value: 6.137e-07
# this is what the raw data look like
plot(x = newDF$deltaLoad, y = newDF$deltaFrq2)

# this package plots the partial residuals
crPlot(m7, variable = "deltaLoad")

# plot with ggplot2
theme_set(theme_classic()) 

# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaArcL, color = deltaLoad)) + 
     geom_point()

ggplot(newDF, aes(x= deltaLoad, y = deltaArcL, color = size_pc1)) + 
     geom_point()

# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) + 
     geom_point() + 
     labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on arcL^2 \n while subtracing affect of bee size") + 
     stat_smooth(method = 'lm', se = FALSE)


Takeaways from model 1:

  1. A higher deltaload (i.e. a relatively larger “high” load) causes bees to have a larger change in arclength squared. Another way to say this is that increasing the change in load by 1 gram causes an increase in the change in arclength^2 by 0.002059 degrees.
  2. The larger the bee, the lower the change in arcLength^2 (association, rather than causation), while holding the change in load constant.
  3. We do not have enough evidence to notice any non-linear relationships in this model.

See how \(\Delta\)load affects \(\Delta\)frequency, accounting for bee size and treatment order

# fit full model
m1 <- lm(deltaFrq2 ~  (deltaLoad +  size_pc1 + order_1)^2  + 
              size_pc1_2  + deltaLoad2, data = newDF)
summary(m1)
## 
## Call:
## lm(formula = deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 + 
##     size_pc1_2 + deltaLoad2, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3981.4 -1255.9  -170.3  1165.6  3552.9 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4195.97     909.39   4.614  0.00015 ***
## deltaLoad                     -106368.80   82390.69  -1.291  0.21073    
## size_pc1                          666.06     883.26   0.754  0.45917    
## order_1loadedSecond             -4376.19    1383.48  -3.163  0.00469 ** 
## size_pc1_2                       -120.74     343.06  -0.352  0.72839    
## deltaLoad2                    -469149.49 3053199.27  -0.154  0.87935    
## deltaLoad:size_pc1               6902.45   58657.40   0.118  0.90744    
## deltaLoad:order_1loadedSecond   62840.01  155625.35   0.404  0.69045    
## size_pc1:order_1loadedSecond       54.39    1467.95   0.037  0.97079    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2119 on 21 degrees of freedom
## Multiple R-squared:  0.5048, Adjusted R-squared:  0.3162 
## F-statistic: 2.676 on 8 and 21 DF,  p-value: 0.03372
car::vif(m1) # some serious multicollinearity
##          deltaLoad           size_pc1            order_1 
##          13.203325          19.003708           3.182227 
##         size_pc1_2         deltaLoad2 deltaLoad:size_pc1 
##          18.740891          10.117446          25.115087 
##  deltaLoad:order_1   size_pc1:order_1 
##          11.384208          23.301468
m2 <- update(m1, .~. - size_pc1:order_1)
anova(m1, m2)
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 + 
##     deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:size_pc1 + deltaLoad:order_1
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     21 94310333                           
## 2     22 94316499 -1   -6165.2 0.0014 0.9708
summary(m2)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + deltaLoad:size_pc1 + deltaLoad:order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3973.1 -1259.6  -176.2  1149.9  3556.7 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4205.4      852.7   4.932 6.21e-05 ***
## deltaLoad                     -108212.1    64166.1  -1.686  0.10584    
## size_pc1                          695.7      367.7   1.892  0.07176 .  
## order_1loadedSecond             -4386.5     1324.0  -3.313  0.00316 ** 
## size_pc1_2                       -109.9      176.0  -0.625  0.53864    
## deltaLoad2                    -417980.4  2660489.9  -0.157  0.87659    
## deltaLoad:size_pc1               5130.1    33170.6   0.155  0.87850    
## deltaLoad:order_1loadedSecond   66238.1   122843.3   0.539  0.59516    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2071 on 22 degrees of freedom
## Multiple R-squared:  0.5048, Adjusted R-squared:  0.3472 
## F-statistic: 3.203 on 7 and 22 DF,  p-value: 0.01701
m3 <- update(m2, .~. - deltaLoad:size_pc1)

anova(m2, m3)
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:size_pc1 + deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     22 94316499                           
## 2     23 94419044 -1   -102545 0.0239 0.8785
summary(m3)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + deltaLoad:order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3910.3 -1244.0  -168.5  1173.4  3503.8 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4190.99     829.34   5.053 4.09e-05 ***
## deltaLoad                     -108478.92   62767.13  -1.728  0.09734 .  
## size_pc1                          702.89     356.93   1.969  0.06108 .  
## order_1loadedSecond             -4434.40    1259.70  -3.520  0.00184 ** 
## size_pc1_2                        -88.67     107.58  -0.824  0.41829    
## deltaLoad2                    -215246.26 2265485.45  -0.095  0.92513    
## deltaLoad:order_1loadedSecond   61643.10  116639.62   0.528  0.60222    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2026 on 23 degrees of freedom
## Multiple R-squared:  0.5042, Adjusted R-squared:  0.3749 
## F-statistic: 3.899 on 6 and 23 DF,  p-value: 0.00785
m4 <- update(m3, .~. - deltaLoad:order_1)
anova(m3, m4) 
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     23 94419044                           
## 2     24 95565635 -1  -1146591 0.2793 0.6022
summary(m4)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3667.1 -1308.7   -52.9  1317.9  3482.6 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          4.330e+03  7.749e+02   5.587 9.48e-06 ***
## deltaLoad           -8.589e+04  4.527e+04  -1.897 0.069889 .  
## size_pc1             7.352e+02  3.463e+02   2.123 0.044282 *  
## order_1loadedSecond -4.794e+03  1.044e+03  -4.594 0.000117 ***
## size_pc1_2          -5.195e+01  8.090e+01  -0.642 0.526837    
## deltaLoad2          -1.250e+06  1.124e+06  -1.112 0.277131    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1995 on 24 degrees of freedom
## Multiple R-squared:  0.4982, Adjusted R-squared:  0.3937 
## F-statistic: 4.766 on 5 and 24 DF,  p-value: 0.00364
m5 <- update(m4, .~. - size_pc1_2)
anova(m4,m5) 
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     24 95565635                           
## 2     25 97207799 -1  -1642165 0.4124 0.5268
summary(m5)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3591.2 -1388.0  -154.6  1391.9  3336.4 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             4166.0      723.2   5.760 5.29e-06 ***
## deltaLoad             -80714.2    44020.7  -1.834 0.078654 .  
## size_pc1                 732.5      342.2   2.141 0.042251 *  
## order_1loadedSecond    -4719.5     1024.9  -4.605 0.000104 ***
## deltaLoad2          -1475818.2  1054446.2  -1.400 0.173916    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1972 on 25 degrees of freedom
## Multiple R-squared:  0.4896, Adjusted R-squared:  0.4079 
## F-statistic: 5.995 on 4 and 25 DF,  p-value: 0.00159
m6 <- update(m5, .~.  - deltaLoad2)
anova(m6, m5) 
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1     26 104824695                           
## 2     25  97207799  1   7616896 1.9589 0.1739
summary(m6)
## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3735.3      666.5   5.605 6.87e-06 ***
## deltaLoad           -105594.7    41007.6  -2.575 0.016065 *  
## size_pc1                861.4      335.6   2.566 0.016382 *  
## order_1loadedSecond   -4717.7     1043.6  -4.521 0.000119 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared:  0.4496, Adjusted R-squared:  0.3861 
## F-statistic:  7.08 on 3 and 26 DF,  p-value: 0.001244
m7 <- update(m6, .~. - size_pc1)
anova(m7, m6) # p-values for size
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1     27 131380648                              
## 2     26 104824695  1  26555953 6.5868 0.01638 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m8 <- update(m6, .~. - deltaLoad)
anova(m6, m8) # p-value for deltaLoad
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ size_pc1 + order_1
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1     26 104824695                              
## 2     27 131557512 -1 -26732816 6.6306 0.01607 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m9 <- update(m6, .~. - order_1)
anova(m6, m9) # p-value for order
## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1     26 104824695                                  
## 2     27 187214161 -1 -82389465 20.435 0.0001192 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit final model with non-centered deltaLoad
m10 <- update(m6, .~. - deltaLoad + dLoad_nonCent)
summary(m10) # final model for paper
## 
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           11857.7     3586.6   3.306 0.002766 ** 
## size_pc1                861.4      335.6   2.566 0.016382 *  
## order_1loadedSecond   -4717.7     1043.6  -4.521 0.000119 ***
## dLoad_nonCent       -105594.7    41007.6  -2.575 0.016065 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared:  0.4496, Adjusted R-squared:  0.3861 
## F-statistic:  7.08 on 3 and 26 DF,  p-value: 0.001244

model diagnostics

par(mfrow = c(2,2))
plot(m10, which = 1:4) # no glaring violations

par(mfrow = c(1,1))

car::vif(m10) # vif is a little high
##      size_pc1       order_1 dLoad_nonCent 
##      3.056456      2.017003      3.643394

model visualization

# calculate partial residuals for deltaLoad
# these are the  residuals, minus the effect of detlaLoad
y <- residuals(m10, type = 'partial')[, "dLoad_nonCent"]

# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)

# double check to make sure the slope for partial residual plots are the 
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent)) 
## 
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             8122       1631   4.980 2.93e-05 ***
## newDF$dLoad_nonCent  -105595      20702  -5.101 2.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1935 on 28 degrees of freedom
## Multiple R-squared:  0.4816, Adjusted R-squared:  0.4631 
## F-statistic: 26.02 on 1 and 28 DF,  p-value: 2.106e-05
# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaFrq2)

# this package plots the partial residuals
crPlot(m10, variable = "dLoad_nonCent")

# plot with ggplot2
# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaFrq2, color = deltaLoad, shape = order_1)) + 
     geom_point()

ggplot(newDF, aes(x= deltaLoad, y = deltaFrq2, color = size_pc1, shape = order_1)) + 
     geom_point()

# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) + 
     geom_point() + 
     labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on freq^2 \n while subtracing affect of bee size and order") + 
     stat_smooth(method = 'lm', se = FALSE)

# partial residuals for order

y <- residuals(m10, type = 'partial')[, "order_1"]

ggplot(newDF, aes(x= order_1, y = y)) + 
     geom_boxplot() + 
     labs(x = "order", y = "partial residuals for order \n i.e. order effect on freq^2 \n while subtracing affect of bee size and load") + 
     stat_smooth(method = 'lm', se = FALSE)

# holding bee size and delta load constant, if a bee was loaded second, (confusing, huh?) then it would have a much lower freq^2 than if it was loaded first

summary(m10)
## 
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           11857.7     3586.6   3.306 0.002766 ** 
## size_pc1                861.4      335.6   2.566 0.016382 *  
## order_1loadedSecond   -4717.7     1043.6  -4.521 0.000119 ***
## dLoad_nonCent       -105594.7    41007.6  -2.575 0.016065 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared:  0.4496, Adjusted R-squared:  0.3861 
## F-statistic:  7.08 on 3 and 26 DF,  p-value: 0.001244

Takeaways from model 2:

  1. Order, size, and deltaLoad are associated with a change in frequency^2. We find no evidence of non-linear relationships or interactions.
  2. Holding size and order constant, we find that a larger deltaLoad causes a decrease in deltaFrequency^2
  3. Holding other variables constant, an increase in bee size is associated with a larger deltaFrequency^2.
  4. Holding other variables constant, if the bee was loaded second, then they had a lower deltaFrequency^2 than if they were loaded first.
# make some predictions
predDF <- data.frame(size_pc1 = 0, order_1 = c("loadedSecond", "loadedFirst"), dLoad_nonCent = mean(newDF$dLoad_nonCent))

predDF$pred_deltaF2 <- predict(m10, predDF)

predDF
##   size_pc1      order_1 dLoad_nonCent pred_deltaF2
## 1        0 loadedSecond       0.07692     -982.314
## 2        0  loadedFirst       0.07692     3735.337

2. We check to see how \(\Delta\)frequency and \(\Delta\)ArcLength are associated with \(\Delta\)MetabolicRate, while accounting for \(\Delta\)load

mm1 <- lm(deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2, data = newDF)
car::vif(mm1)
##  deltaFrq2 deltaArcL2  deltaLoad   size_pc1 deltaLoad2 
##   1.060059   1.852143   3.844322   2.679591   1.495113
summary(mm1)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + 
##     size_pc1 + deltaLoad2, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.11198 -0.64060 -0.01055  0.42457  2.97948 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.941e+00  6.182e-01   3.140  0.00444 ** 
## deltaFrq2    3.543e-04  6.888e-05   5.144 2.89e-05 ***
## deltaArcL2  -2.573e+03  5.900e+03  -0.436  0.66659    
## deltaLoad    2.368e+01  1.937e+01   1.223  0.23337    
## size_pc1    -5.482e-02  1.445e-01  -0.379  0.70776    
## deltaLoad2   4.440e+02  5.114e+02   0.868  0.39381    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9233 on 24 degrees of freedom
## Multiple R-squared:  0.5937, Adjusted R-squared:  0.5091 
## F-statistic: 7.014 on 5 and 24 DF,  p-value: 0.0003629
mm2 <- update(mm1, .~. - size_pc1)
anova(mm1, mm2)
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     24 20.460                           
## 2     25 20.582 -1  -0.12268 0.1439 0.7078
summary(mm2)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + 
##     deltaLoad2, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06884 -0.60009 -0.03648  0.36122  3.09536 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.821e+00  5.223e-01   3.487  0.00182 ** 
## deltaFrq2    3.571e-04  6.730e-05   5.306 1.69e-05 ***
## deltaArcL2  -1.543e+03  5.147e+03  -0.300  0.76682    
## deltaLoad    1.788e+01  1.168e+01   1.530  0.13849    
## deltaLoad2   4.891e+02  4.888e+02   1.001  0.32665    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9074 on 25 degrees of freedom
## Multiple R-squared:  0.5913, Adjusted R-squared:  0.5259 
## F-statistic: 9.041 on 4 and 25 DF,  p-value: 0.0001166
mm3 <- update(mm2, .~. - deltaArcL2)
anova(mm3, mm2)
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     26 20.656                           
## 2     25 20.582  1   0.07399 0.0899 0.7668
summary(mm3)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2, 
##     data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.05850 -0.58135 -0.00991  0.37733  3.10956 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.682e+00  2.327e-01   7.228 1.12e-07 ***
## deltaFrq2   3.560e-04  6.602e-05   5.393 1.19e-05 ***
## deltaLoad   1.667e+01  1.077e+01   1.548    0.134    
## deltaLoad2  4.422e+02  4.550e+02   0.972    0.340    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8913 on 26 degrees of freedom
## Multiple R-squared:  0.5898, Adjusted R-squared:  0.5425 
## F-statistic: 12.46 on 3 and 26 DF,  p-value: 3.065e-05
mm4 <- update(mm3, .~. - deltaLoad2)
anova(mm3, mm4)
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     26 20.656                           
## 2     27 21.407 -1  -0.75053 0.9447   0.34
summary(mm4)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.824e+00  1.808e-01  10.086 1.18e-10 ***
## deltaFrq2   3.451e-04  6.498e-05   5.310 1.32e-05 ***
## deltaLoad   2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06
mm5 <- update(mm4, .~. - deltaLoad) 
anova(mm5, mm4) # p-value for load
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     28 25.352                              
## 2     27 21.407  1    3.9448 4.9754 0.03421 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mm6 <- update(mm4, .~. - deltaFrq2)
anova(mm4, mm6) # p-value for deltafrq2
## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad
## Model 2: deltaMetR ~ deltaLoad
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     27 21.407                                  
## 2     28 43.766 -1   -22.359 28.201 1.324e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit model with non-centered load
mm7 <- update(mm4, .~. - deltaLoad + dLoad_nonCent)

summary(mm7) # final model for paper
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.777e-01  7.507e-01   0.237   0.8146    
## deltaFrq2     3.451e-04  6.498e-05   5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06

model diagnostics

car::vif(mm7)
##     deltaFrq2 dLoad_nonCent 
##      1.014369      1.014369
par(mfrow = c(2,2))
plot(mm7, which = 1:4) # no glaring violations

par(mfrow = c(1,1))

model visualization

# calculate partial residuals for deltaLoad
# these are the  residuals, minus the effect of detlaLoad
y <- residuals(mm7, type = 'partial')[, "dLoad_nonCent"]

# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)

# double check to make sure the slope for partial residual plots are the 
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent)) 
## 
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)  
## (Intercept)          -1.6463     0.7371  -2.233   0.0337 *
## newDF$dLoad_nonCent  21.4030     9.3554   2.288   0.0299 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8744 on 28 degrees of freedom
## Multiple R-squared:  0.1575, Adjusted R-squared:  0.1274 
## F-statistic: 5.234 on 1 and 28 DF,  p-value: 0.02991
summary(mm7)
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.777e-01  7.507e-01   0.237   0.8146    
## deltaFrq2     3.451e-04  6.498e-05   5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06
# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaMetR)

# this package plots the partial residuals
crPlot(mm7, variable = "dLoad_nonCent")

crPlot(mm7, variable = "deltaFrq2")

summary(mm7) # final model for paper
## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.777e-01  7.507e-01   0.237   0.8146    
## deltaFrq2     3.451e-04  6.498e-05   5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06

Takeaways for model 3:

  1. Deltafreq^2 and deltaLoad are both associated with deltaMetR.
  2. Holding deltafreq constant (not necessarily at 0), an increase in deltaload is associated with an increase in deltametabolicRate
  3. Holding deltaLoad constant (not holding load constant, and again, not holding deltaload necessarily at 0), an increase in deltaFreq^2 is associated with an increase in deltaMetabolic rate.
  4. We found no evidence to suggest that a change in deltaArcL^2 is associated with deltametabolicRate. (this is not saying that arcLength doesn’t affect metabolic rate).

Refref: Do regression for Average percent loading vs. delta metabolic rate / 1% load etc. It’s the last page of the document Susie sent (3 different response variables).

Covariates: avg % load, order, bee size.

(DeltaMetRate / avgPercLoad) ~ avgPercLoad + order_1 + size_pc1

newDF <- within(newDF, {dmr_apl = deltaMetR / avgPercLoad})

mod1 <- lm(dmr_apl  ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1)
## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.029235 -0.014142 -0.003264  0.008463  0.053034 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.1404689  0.0176079   7.978 1.86e-08 ***
## avgPercLoad         -0.0015015  0.0003086  -4.865 4.80e-05 ***
## order_1loadedSecond -0.0212737  0.0073765  -2.884  0.00779 ** 
## size_pc1             0.0012171  0.0019858   0.613  0.54527    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01965 on 26 degrees of freedom
## Multiple R-squared:  0.6107, Adjusted R-squared:  0.5657 
## F-statistic: 13.59 on 3 and 26 DF,  p-value: 1.58e-05
mod2 <- update(mod1, .~. - size_pc1)
anova(mod1, mod2) # remove size_pc1
## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad + order_1 + size_pc1
## Model 2: dmr_apl ~ avgPercLoad + order_1
##   Res.Df      RSS Df   Sum of Sq      F Pr(>F)
## 1     26 0.010044                             
## 2     27 0.010189 -1 -0.00014511 0.3756 0.5453
summary(mod2) # final model for paper
## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad + order_1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030989 -0.013911 -0.001532  0.009837  0.050990 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.1431985  0.0168372   8.505 4.06e-09 ***
## avgPercLoad         -0.0015573  0.0002914  -5.344 1.21e-05 ***
## order_1loadedSecond -0.0204551  0.0071702  -2.853  0.00822 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01943 on 27 degrees of freedom
## Multiple R-squared:  0.605,  Adjusted R-squared:  0.5758 
## F-statistic: 20.68 on 2 and 27 DF,  p-value: 3.578e-06
# get p-values for avgPercLoad and order
anova(mod2,  update(mod2, .~. - order_1)) # p-value for order
## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad + order_1
## Model 2: dmr_apl ~ avgPercLoad
##   Res.Df      RSS Df  Sum of Sq      F   Pr(>F)   
## 1     27 0.010189                                 
## 2     28 0.013260 -1 -0.0030713 8.1384 0.008216 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod2, update(mod2, .~. - avgPercLoad)) # p-value for avg perc load
## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad + order_1
## Model 2: dmr_apl ~ order_1
##   Res.Df      RSS Df Sum of Sq      F    Pr(>F)    
## 1     27 0.010189                                  
## 2     28 0.020965 -1 -0.010776 28.554 1.211e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# plot(I(deltaMetR)  ~ avgPercLoad, data = newDF)
# plot(I(deltaMetR / avgPercLoad)  ~ deltaMetR, data = newDF)

car::scatterplotMatrix(newDF[, c("avgPercLoad", "deltaMetR", "size_pc1")])

car::vif(mod1)  # looks fine
## avgPercLoad     order_1    size_pc1 
##    1.114454    1.051699    1.116654
car::vif(mod2)
## avgPercLoad     order_1 
##    1.017216    1.017216
par(mfrow = c(2,2))
plot(mod2) # possible nonlinear trend in residuals

plot(mod1, which = 4)

# update model to add a non-linear term
newDF <- within(newDF, {avgPercLoad_cent = as.numeric(scale(newDF$avgPercLoad, center = TRUE, scale = FALSE))})
newDF$apl_cent2 <- with(newDF, avgPercLoad_cent^2)

m11 <- lm(dmr_apl  ~ avgPercLoad_cent  + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m11)
## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + 
##     size_pc1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030273 -0.010355 -0.003798  0.006817  0.055061 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          5.247e-02  5.796e-03   9.053 2.29e-09 ***
## avgPercLoad_cent    -1.553e-03  3.088e-04  -5.030 3.45e-05 ***
## apl_cent2            2.772e-05  2.280e-05   1.216  0.23544    
## order_1loadedSecond -2.364e-02  7.564e-03  -3.125  0.00446 ** 
## size_pc1             1.321e-03  1.970e-03   0.671  0.50864    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01948 on 25 degrees of freedom
## Multiple R-squared:  0.6324, Adjusted R-squared:  0.5736 
## F-statistic: 10.75 on 4 and 25 DF,  p-value: 3.289e-05
car::vif(m11) # vif is much better with centered variables
## avgPercLoad_cent        apl_cent2          order_1         size_pc1 
##         1.135857         1.108436         1.126036         1.118753
m11a <- update(m11, .~. - size_pc1)
anova(m11, m11a) # remove size_pc1
## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1
##   Res.Df       RSS Df   Sum of Sq      F Pr(>F)
## 1     25 0.0094835                             
## 2     26 0.0096541 -1 -0.00017058 0.4497 0.5086
m11b <- update(m11a, .~. - apl_cent2)
anova(m11a, m11b) # drop squared term
## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1
## Model 2: dmr_apl ~ avgPercLoad_cent + order_1
##   Res.Df       RSS Df   Sum of Sq      F Pr(>F)
## 1     26 0.0096541                             
## 2     27 0.0101893 -1 -0.00053522 1.4414 0.2407
summary(m11b)
## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad_cent + order_1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030989 -0.013911 -0.001532  0.009837  0.050990 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.0549528  0.0052157  10.536 4.58e-11 ***
## avgPercLoad_cent    -0.0015573  0.0002914  -5.344 1.21e-05 ***
## order_1loadedSecond -0.0204551  0.0071702  -2.853  0.00822 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01943 on 27 degrees of freedom
## Multiple R-squared:  0.605,  Adjusted R-squared:  0.5758 
## F-statistic: 20.68 on 2 and 27 DF,  p-value: 3.578e-06
par(mfrow = c(2,2))

plot(m11b) # residuals look much better

plot(m11b, which = 4)

summary(m11b)
## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad_cent + order_1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030989 -0.013911 -0.001532  0.009837  0.050990 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.0549528  0.0052157  10.536 4.58e-11 ***
## avgPercLoad_cent    -0.0015573  0.0002914  -5.344 1.21e-05 ***
## order_1loadedSecond -0.0204551  0.0071702  -2.853  0.00822 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01943 on 27 degrees of freedom
## Multiple R-squared:  0.605,  Adjusted R-squared:  0.5758 
## F-statistic: 20.68 on 2 and 27 DF,  p-value: 3.578e-06
## visualize model for deltametRate/avgPercLoading

par(mfrow = c(1,1))

plot(dmr_apl  ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)

ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200), 
                  order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))

ndf$apl_cent2 <- ndf$avgPercLoad_cent^2

ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)
predDF <- data.frame(preds = predict(mod2, newdata = ndf), ndf)

ggplot(newDF, aes(x = avgPercLoad, y = dmr_apl)) + 
     geom_point(aes(color = order_1), , shape = 17) + 
     geom_line(data = predDF, aes(x = avgPercLoad, y = preds, color = order_1)) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in Metabolic Rate (mL CO2 / hr) / \n Average Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8))

ggsave("dmr_apl.pdf", width = 5, height = 4)

(DeltaFrq^2/avgPercload) ~ avgPercLoad + order_1 + size_pc1

newDF <- within(newDF, {df2_apl = deltaFrq2 / avgPercLoad})
par(mfrow = c(1,1))
plot(df2_apl ~ avgPercLoad, newDF)

mod1_f <- lm(df2_apl  ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_f)
## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.153 -25.697   3.939  22.568  51.988 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         221.4973    29.3034   7.559 5.04e-08 ***
## avgPercLoad          -3.0335     0.5136  -5.906 3.14e-06 ***
## order_1loadedSecond -37.1107    12.2762  -3.023  0.00557 ** 
## size_pc1             -2.6032     3.3048  -0.788  0.43800    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.71 on 26 degrees of freedom
## Multiple R-squared:  0.667,  Adjusted R-squared:  0.6286 
## F-statistic: 17.36 on 3 and 26 DF,  p-value: 2.15e-06
mod2_f <- update(mod1_f, .~. - size_pc1)
anova(mod1_f, mod2_f) # remove size_pc1
## Analysis of Variance Table
## 
## Model 1: df2_apl ~ avgPercLoad + order_1 + size_pc1
## Model 2: df2_apl ~ avgPercLoad + order_1
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     26 27819                           
## 2     27 28482 -1   -663.86 0.6205  0.438
summary(mod2_f) # not final model for paper
## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad + order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -54.479 -27.491   1.816  21.916  49.750 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         215.6591    28.1506   7.661 3.07e-08 ***
## avgPercLoad          -2.9140     0.4873  -5.980 2.23e-06 ***
## order_1loadedSecond -38.8616    11.9881  -3.242  0.00315 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.48 on 27 degrees of freedom
## Multiple R-squared:  0.6591, Adjusted R-squared:  0.6338 
## F-statistic:  26.1 on 2 and 27 DF,  p-value: 4.904e-07
# get p-values for avgPercLoad 
anova(mod2_f,  update(mod2_f, .~. - order_1)) # p-value for order
## Analysis of Variance Table
## 
## Model 1: df2_apl ~ avgPercLoad + order_1
## Model 2: df2_apl ~ avgPercLoad
##   Res.Df   RSS Df Sum of Sq      F   Pr(>F)   
## 1     27 28482                                
## 2     28 39568 -1    -11086 10.508 0.003152 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow = c(2,2))
plot(mod2_f) # non-linearity in residuals

plot(mod2_f, which = 4)


# update model to add a non-linear term
m22 <- lm(df2_apl  ~ avgPercLoad_cent  + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m22)
## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + 
##     size_pc1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -35.430 -16.698  -1.678  19.332  43.880 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          35.70197    7.30757   4.886 5.01e-05 ***
## avgPercLoad_cent     -3.27920    0.38928  -8.424 9.05e-09 ***
## apl_cent2             0.13215    0.02874   4.597 0.000106 ***
## order_1loadedSecond -48.37471    9.53594  -5.073 3.09e-05 ***
## size_pc1             -2.10870    2.48324  -0.849 0.403845    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.56 on 25 degrees of freedom
## Multiple R-squared:  0.8196, Adjusted R-squared:  0.7907 
## F-statistic: 28.39 on 4 and 25 DF,  p-value: 5.684e-09
car::vif(m22) # vif is much better with centered variables
## avgPercLoad_cent        apl_cent2          order_1         size_pc1 
##         1.135857         1.108436         1.126036         1.118753
m22a <- update(m22, .~. - size_pc1)
anova(m22, m22a) # remove size_pc1
## Analysis of Variance Table
## 
## Model 1: df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     25 15074                           
## 2     26 15509 -1    -434.8 0.7211 0.4038
summary(m22a)
## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.409 -16.168  -2.024  16.653  47.040 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          36.34581    7.22905   5.028 3.12e-05 ***
## avgPercLoad_cent     -3.18454    0.37097  -8.584 4.60e-09 ***
## apl_cent2             0.13321    0.02856   4.664 8.17e-05 ***
## order_1loadedSecond -49.88054    9.31922  -5.352 1.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.42 on 26 degrees of freedom
## Multiple R-squared:  0.8144, Adjusted R-squared:  0.793 
## F-statistic: 38.02 on 3 and 26 DF,  p-value: 1.183e-09
par(mfrow = c(2,2))

plot(m22a) # residuals look much better

plot(m22a, which = 4)

summary(m22a) # final model for paper (though we could also report the non-centered values)
## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.409 -16.168  -2.024  16.653  47.040 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          36.34581    7.22905   5.028 3.12e-05 ***
## avgPercLoad_cent     -3.18454    0.37097  -8.584 4.60e-09 ***
## apl_cent2             0.13321    0.02856   4.664 8.17e-05 ***
## order_1loadedSecond -49.88054    9.31922  -5.352 1.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.42 on 26 degrees of freedom
## Multiple R-squared:  0.8144, Adjusted R-squared:  0.793 
## F-statistic: 38.02 on 3 and 26 DF,  p-value: 1.183e-09
## visualize model for deltametRate/avgPercLoading

ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200), 
                  order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))

ndf$apl_cent2 <- ndf$avgPercLoad_cent^2

ndf$apl <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)

predDF <- data.frame(preds = predict(m22a, newdata = ndf, se = TRUE), ndf)

par(mfrow = c(1,1))

plot(df2_apl  ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)

ggplot(newDF, aes(x = avgPercLoad, y = df2_apl)) + 
     geom_point(aes(color = order_1), , shape = 17) + 
     geom_line(data = predDF, aes(x = apl, y = preds.fit, color = order_1)) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in wingbeat freq^2 (hz^2) / \n Average Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8)) 

ggsave("df2_apl.pdf", width = 5, height = 4)

(DeltaArcL^2/avgPercload) ~ avgPercLoad + order_1 + size_pc1

newDF <- within(newDF, {da2_apl = deltaArcL2 / avgPercLoad})

par(mfrow = c(1,1))
plot(da2_apl ~ avgPercLoad, ylab = c("deltaArcL^2 / avgPercLoad"), newDF)

mod1_a <- lm(da2_apl  ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_a)
## 
## Call:
## lm(formula = da2_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.023e-06 -5.865e-07  7.780e-08  4.151e-07  1.346e-06 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          3.321e-06  6.130e-07   5.418 1.12e-05 ***
## avgPercLoad         -2.018e-08  1.074e-08  -1.878   0.0717 .  
## order_1loadedSecond -6.484e-07  2.568e-07  -2.525   0.0180 *  
## size_pc1             1.764e-08  6.913e-08   0.255   0.8006    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.843e-07 on 26 degrees of freedom
## Multiple R-squared:  0.3183, Adjusted R-squared:  0.2397 
## F-statistic: 4.047 on 3 and 26 DF,  p-value: 0.01737
mod2_a <- update(mod1_a, .~. - size_pc1)
anova(mod1_a, mod2_a) # remove size_pc1
## Analysis of Variance Table
## 
## Model 1: da2_apl ~ avgPercLoad + order_1 + size_pc1
## Model 2: da2_apl ~ avgPercLoad + order_1
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     26 1.2174e-11                             
## 2     27 1.2204e-11 -1 -3.0501e-14 0.0651 0.8006
summary(mod2_a) # final model for paper
## 
## Call:
## lm(formula = da2_apl ~ avgPercLoad + order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.033e-06 -5.437e-07  9.477e-08  4.213e-07  1.332e-06 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          3.361e-06  5.827e-07   5.768 3.91e-06 ***
## avgPercLoad         -2.099e-08  1.009e-08  -2.081   0.0471 *  
## order_1loadedSecond -6.365e-07  2.482e-07  -2.565   0.0162 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.723e-07 on 27 degrees of freedom
## Multiple R-squared:  0.3166, Adjusted R-squared:  0.266 
## F-statistic: 6.255 on 2 and 27 DF,  p-value: 0.005861
# get p-values for avgPercLoad and order
anova(mod2_a,  update(mod2_a, .~. - order_1)) # p-value for order
## Analysis of Variance Table
## 
## Model 1: da2_apl ~ avgPercLoad + order_1
## Model 2: da2_apl ~ avgPercLoad
##   Res.Df        RSS Df   Sum of Sq    F  Pr(>F)  
## 1     27 1.2204e-11                              
## 2     28 1.5178e-11 -1 -2.9742e-12 6.58 0.01619 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod2_a, update(mod2_a, .~. - avgPercLoad)) # p-value for avg perc load
## Analysis of Variance Table
## 
## Model 1: da2_apl ~ avgPercLoad + order_1
## Model 2: da2_apl ~ order_1
##   Res.Df        RSS Df   Sum of Sq     F  Pr(>F)  
## 1     27 1.2204e-11                               
## 2     28 1.4161e-11 -1 -1.9567e-12 4.329 0.04708 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow = c(2,2))
plot(mod2_a)

plot(mod2_a, which = 4)

# visualize model

ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200), 
                  order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))

ndf$apl_cent2 <- ndf$avgPercLoad_cent^2

ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)

predDF <- data.frame(preds = predict(mod2_a, newdata = ndf, se = TRUE), ndf)

par(mfrow = c(1,1))

plot(da2_apl  ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)

ggplot(newDF, aes(x = avgPercLoad, y = da2_apl)) + 
     geom_point(aes(color = order_1), shape = 17) + 
     geom_line(data = predDF, aes(x = avgPercLoad, y = preds.fit, color = order_1)) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in arc length^2 (radians^2) / \n Average Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8)) 

ggsave("da2_apl.pdf", width = 5, height = 4)

Session Info

sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] plyr_1.8.4         viridis_0.3.4      visreg_2.3-0      
##  [4] effects_3.1-2      sjPlot_2.2.1       influence.ME_0.9-8
##  [7] reshape2_1.4.2     MASS_7.3-45        gsheet_0.4.2      
## [10] lme4_1.1-12        Matrix_1.2-8       car_2.1-4         
## [13] ggplot2_2.2.1     
## 
## loaded via a namespace (and not attached):
##  [1] tidyr_0.6.1        splines_3.3.2      merTools_0.3.0    
##  [4] modelr_0.1.0       shiny_1.0.0        assertthat_0.1    
##  [7] stats4_3.3.2       coin_1.1-3         yaml_2.1.14       
## [10] backports_1.0.5    lattice_0.20-34    quantreg_5.29     
## [13] digest_0.6.12      minqa_1.2.4        colorspace_1.3-2  
## [16] sandwich_2.3-4     htmltools_0.3.5    httpuv_1.3.3      
## [19] psych_1.6.12       broom_0.4.1        SparseM_1.74      
## [22] haven_1.0.0        purrr_0.2.2        xtable_1.8-2      
## [25] mvtnorm_1.0-5      scales_0.4.1       stringdist_0.9.4.4
## [28] MatrixModels_0.4-1 arm_1.9-3          tibble_1.2        
## [31] mgcv_1.8-16        DT_0.2             TH.data_1.0-8     
## [34] nnet_7.3-12        lazyeval_0.2.0     pbkrtest_0.4-6    
## [37] mnormt_1.5-5       survival_2.40-1    magrittr_1.5      
## [40] mime_0.5           evaluate_0.10      nlme_3.1-130      
## [43] foreign_0.8-67     tools_3.3.2        multcomp_1.4-6    
## [46] stringr_1.1.0      munsell_0.4.3      blme_1.0-4        
## [49] grid_3.3.2         nloptr_1.0.4       htmlwidgets_0.8   
## [52] labeling_0.3       rmarkdown_1.3      gtable_0.2.0      
## [55] codetools_0.2-15   curl_2.3           abind_1.4-5       
## [58] sjstats_0.7.1      DBI_0.5-1          sjmisc_2.2.1      
## [61] R6_2.2.0           gridExtra_2.2.1    zoo_1.7-14        
## [64] knitr_1.15.1       dplyr_0.5.0        rprojroot_1.2     
## [67] readr_1.0.0        modeltools_0.2-21  stringi_1.1.2     
## [70] parallel_3.3.2     Rcpp_0.12.9        coda_0.19-1       
## [73] lmtest_0.9-34
Sys.time()
## [1] "2017-03-16 10:12:29 EDT"